/*
 * beamMacro.cc
 *
 *  Created on: Mar 14, 2011
 *      Author: dblyth
 */

#include "iostream"
#include "cmath"
using namespace std;

int main(int argc, char** argv) {

	const double pi = 3.14159265;
	const double w = 127;
	const double L = 152.4;
	const double dAngle = pi / 16.0;
	const double dLength = 20;

	double lowerSpaceH = 161;
	double aerogelSpaceH = 112;
	double mainBoxThickness = 1.5;

	double h = 312;

	double minX = 0;
	double maxX = L;
	double x = minX;

	for (int ix = 0; x <= maxX; ix++) {
		double minY = 0;
		double maxY = w;
		double y = minY;

		for (int iy = 0; y <= maxY; iy++) {
			double minPhi = atan((w - y) / (L - x));
			double maxPhi = pi + atan((w - y) / (-x));
			double phi = minPhi;

			for (int iphi = 0; phi <= maxPhi; iphi++) {
				double minTheta = 0;
				double maxTheta = atan((w - y) / cos(phi) / h / tan(phi));
				double theta = minTheta;

				for (int itheta = 0; theta <= maxTheta; itheta++) {
					cout << "/gun/position " << y - w / 2 << " " << x - L / 2 << " " << (lowerSpaceH + aerogelSpaceH) / 2 + mainBoxThickness
							<< " mm\n";
					cout << "/gun/direction " << sin(theta) * sin(phi) << " " << sin(theta) * cos(phi) << " " << -cos(theta) << endl;
					cout << "/run/beamOn 1\n";

					theta = itheta * dAngle + minTheta;
				}

				phi = iphi * dAngle + minPhi;
			}

			y = iy * dLength + minY;
		}

		x = ix * dLength + minX;
	}

	x = minX;

	for (int ix = 0; x <= maxX; ix++) {
		double minY = 0;
		double maxY = w;
		double y = minY;

		for (int iy = 0; y <= maxY; iy++) {
			double minPhi = atan((-y) / (L - x));
			double maxPhi = atan((w - y) / (L - x));
			double phi = minPhi;

			for (int iphi = 0; phi <= maxPhi; iphi++) {
				double minTheta = 0;
				double maxTheta = atan((L - x) / cos(phi) / h);
				double theta = minTheta;

				for (int itheta = 0; theta <= maxTheta; itheta++) {
					cout << "/gun/position " << y - w / 2 << " " << x - L / 2 << " " << (lowerSpaceH + aerogelSpaceH) / 2 + mainBoxThickness
							<< " mm\n";
					cout << "/gun/direction " << sin(theta) * sin(phi) << " " << sin(theta) * cos(phi) << " " << -cos(theta) << endl;
					cout << "/run/beamOn 1\n";

					theta = itheta * dAngle + minTheta;
				}

				phi = iphi * dAngle + minPhi;
			}

			y = iy * dLength + minY;
		}

		x = ix * dLength + minX;
	}

	x = minX;

	for (int ix = 0; x <= maxX; ix++) {
		double minY = 0;
		double maxY = w;
		double y = minY;

		for (int iy = 0; y <= maxY; iy++) {
			double minPhi = atan((w - y) / (L - x));
			double maxPhi = pi + atan((w - y) / (-x));
			double phi = minPhi;

			for (int iphi = 0; phi <= maxPhi; iphi++) {
				double minTheta = 0;
				double maxTheta = atan((w - y) / cos(phi) / h / tan(phi));
				double theta = minTheta;

				for (int itheta = 0; theta <= maxTheta; itheta++) {
					cout << "/gun/position " << -y + w / 2 << " " << -x + L / 2 << " " << (lowerSpaceH + aerogelSpaceH) / 2 + mainBoxThickness
							<< " mm\n";
					cout << "/gun/direction " << -sin(theta) * sin(phi) << " " << -sin(theta) * cos(phi) << " " << -cos(theta) << endl;
					cout << "/run/beamOn 1\n";

					theta = itheta * dAngle + minTheta;
				}

				phi = iphi * dAngle + minPhi;
			}

			y = iy * dLength + minY;
		}

		x = ix * dLength + minX;
	}

	x = minX;

	for (int ix = 0; x <= maxX; ix++) {
		double minY = 0;
		double maxY = w;
		double y = minY;

		for (int iy = 0; y <= maxY; iy++) {
			double minPhi = atan((-y) / (L - x));
			double maxPhi = atan((w - y) / (L - x));
			double phi = minPhi;

			for (int iphi = 0; phi <= maxPhi; iphi++) {
				double minTheta = 0;
				double maxTheta = atan((L - x) / cos(phi) / h);
				double theta = minTheta;

				for (int itheta = 0; theta <= maxTheta; itheta++) {
					cout << "/gun/position " << -y + w / 2 << " " << -x + L / 2 << " " << (lowerSpaceH + aerogelSpaceH) / 2 + mainBoxThickness
							<< " mm\n";
					cout << "/gun/direction " << -sin(theta) * sin(phi) << " " << -sin(theta) * cos(phi) << " " << -cos(theta) << endl;
					cout << "/run/beamOn 1\n";

					theta = itheta * dAngle + minTheta;
				}

				phi = iphi * dAngle + minPhi;
			}

			y = iy * dLength + minY;
		}

		x = ix * dLength + minX;
	}

}
